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Remarques sur la stabilite des PMLs Cartesiennes dans les coins 



Resume : Ce travail est une contribution a la comprehension de la question de la stabilite des couches 
absorbantes parfaitement adaptees (PMLs) dans les coins, aux niveaux continu et discret. Des resultats 
de stabiUte sont d'abord presentes pour des PMLs Cartesiennes associ^es a un systeme hyperboUque de 
premier ordre general. Puis, dans le cadre de la formulation du premier ordre pression-vitesse de 1' equation 
des ondes acoustiques, une formulation non splittee des PMLs est discretisee par des elements finis mixtes 
spectraux en espace et des differences finies en temps. On montre, a travers 1' analyse de stabilite de deux 
schemas, conmient un mauvais choix pour la discretisation en temps peut deteriorer la condition de stability 
CFL. Ces resultats de stabiUte sont iUustres par des experiences numeriques. 

Mots-cles : PML, couches absorbantes parfaitement adaptees, stabilite, elements finis, differences finies, 
condition CFL 
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1 Introduction 



The Perfectly Matched Layers (PML) technique, introduced in 1994 by Berenger 111211 for electromagnetic 
problems, is considered as an efficient tool to simulate numerically wave propagation problems stated in 
unbounded domains. In fact, during the last decade, this technique has been intensively applied in a wide 
range of areas (acoustic and electromagnetic problems [15, 16, 29, 30. 33. 36.] . elastodynamics I6..7ll9l [l7ll . 
aeroacoustics among others). 

The mathematical analysis of the well-posedness and stability of the continuous PML models (both in 
the original split formulation and in unsplit formulations) have been addressed in several works. Let us 
cite e.g. m 12, 27 , 3ll for the well-posedness analysis and 

sea 

for the stability analysis. It is now 

well-known that the original split PML model, as well as unsplit PML models, are (at least weakly) stable 
if the original physical model is isotropic. However, if the physical model is anisotropic then the PML 
technique can lead to unstable behaviors (see ^ for more details). 

On the other hand, a natural question is the analysis of the stability of fully discrete schemes used for 
discretizing PMLs. In 1 10], this question has been considered for the discretization of isotropic Maxwell's 
equations. Precisely, it is shown that the Yee scheme, applied for discretizing both the split and the unsplit 
PMLs, is stable under the standard CFL stability condition. This result is obtained for a layer in only one 
direction. When considering Cartesian PMLs with layers in two (or three) directions, it is then natural to 
address the question of stability in comers. It is well known that this question is delicate for Absorbing 



Boundary Conditions (ABCs) and has been considered by several authors (e.g., H19L 13211 . and 112211 where 
long time instabilities for high order ABCs are mentioned). 

The main goal of the present paper is to study if instabilities could be generated from the corner PML 
domains in 2D, as it is the case for some ABCs in corners. The first result is that, on the continuous 
level, PML comers are always stable (even in anisotropic models assuming that the PML parameters for 
each direction are constant and equal). On the discrete level, it is shown that instabilities are related to an 
inadequate time discretization of the auxiliary differential equations in the Cartesian PML formulation. In 
fact, in the context of the isotropic acoustic model, we analyze two different time discretizations whose 
CFL stability conditions are different only in the comer domains. Although both discrete scheme are 
consistent with the continuous model, only one has a CFL condition independent of the PML parameters, 
which coincides with the standard CFL of the scheme in the physical domain. 

Following is the outline of the paper In Section |2] we describe the split Berenger's PML formulation, 
written in Cartesian coordinates for a general first order hyperbolic system in two dimensions. We focus 
our attention on the equation system stated only in a corner domain. It is shown that, in contrast to a layer in 
only one direction, PMLs are always stable in a corner, at least for a constant damping factor In particular, 
it means that, even for anisotropic models, the split PMLs are stable at the continuous level. 

Section|3]is devoted to the introduction of a model problem, the two dimensional wave equation, written 
as a first-order pressure-velocity system. The PML model considered here is the Zhao-Cangellaris unsplit 
formulation. Again, if the PML coefficients are assumed constant in the corner domain (but possibly 
different), it is shown, via energy estimates, that the continuous model is stable. 

The spatial semi-discretization, presented in Section |4] is done with a spectral mixed finite element 
method based on a quadrilateral mesh, used in [2Q]. 

Section |5] is devoted to the time discretization, using explicit second-order finite differences. We first 
consider the scheme used in i20ll . A stability analysis of the scheme shows that the CFL condition is dete- 
riorated in the PML corner, when compared to the CFL condition inside the fluid domain, the deterioration 
depending (in particular) on the damping factor In order to avoid this, a new discretization in time in the 
PML comer is proposed, for which it is shown that the CFL condition remains the same as in the fluid 
domain. These results are illustrated in Section |6]with numerical simulations. 



2 Stability of Berenger's splitted PMLs in corners for a general hy- 
perbolic system 

Notation. Through the rest of the paper, standard notations about functional Sobolev spaces are used 
without explicit definitions, || -11^2 and (•, ■)i,2 denote respectively the L^-norm and the L^-inner product. 
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In UtIi . the authors have shown how to design a spUtted PML model, for a general first-order hyperbolic 
system. This construction has been applied in |9] for analyzing the well-posedness and the stability of 
a layer in one direction. In this section, we focus on the PML equations written in a corner domain. 
Following jgi flTll . we consider the first-order hyperboUc general Cauchy problem in the two dimensional 
free space: 

dtU - A.,d:,U ~ AydyU ^ 

f/(.,0) = C/", 

where J7 is a to dimensional real-valued vector function U : (x, y, t) £ x R+ U (a;, y, t) £ R™, 
and Ay are m x m real-valued symmetric matrices, and the initial data is a m dimensional real-valued 
vector function [/" : {x, y) G R^ i— > U^{x, y) G R™. It is then classical to show that the solution satisfies 
the energy conservation 

4llc/lli.=o. 



dt 



We now introduce the splitted Berenger's PMLs equations in Cartesian coordinates r ill2lll7ll ') using the 
splitting trick for the corner domain: 

Find (U^ , Ijy) solution of the first-order hyperbolic system 

dtU"" + a^U'= - A,d.,U ^0, (1) 

dtUy + GyUy ~ AydyU = 0, (2) 

[/ c/^ + uy, (3) 
f/^(., 0) = {u'-'f, uy{., 0) = {uyf. (4) 

It is straightforward to verify the following 
Theorem 1. Ifdx and Uy are constant and equal to a then the solution {U^ , U^) of system satisfies 

j^\\U\\l.^-2a\\U\\l..<Q. (5) 

In this case, the PML corner is strongly stable in the sense that the solution can be bounded by the initial 
data, uniformly in time: 

3O0, ||C/(i)||^2 < C||C/°||^, , Vt>0. (6) 

Proof. Since the absorbing functions are equal, adding equations ([T]l and (|2]l, and using (O, we see that U 
satisfies the following equation: 

dtU + nU - A.,d:,U - AydyU = Q. (7) 

In this case, the PML comer appears as a "classical" dissipation term. The energy identity (|5]l and the 
estimate (|6]l follow then easily. □ 

The strong stability shows in particular that the solution can not grow linearly in time, therefore there 
is no possible long-time instability coming from the corner, when the damping factor is the same in the two 
directions (contrarely to some high order ABCs |l22|). 

The proof of an analogous energy estimate, for non-constant and different absorbing functions and 
CTy, remains open (see, for instance [ 10]). 

Remark 1. For several models as, for instance. Maxwell's equations or the acoustic scalar wave equation, 
several authors iTil 75 \34il derived the Berenger's PML model in terms of a complex coordinate 



stretching in the frequency domain. In the time domain, this point of view leads to the construction of the 
so-called unsplit PMLs, which do not need the splitting of the unknowns fields, but the introduction of new 
additional unknowns. The main advantage of this formulation is that it preserves the spatial differential 
operators and consequently the original spatial discretizations developed for the physical models. The 
original Berenger's PML can then be reformulated in several forms (e.g. U^ .IL .2L,35,1 ]. depending on the 
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choice of the auxiliary unknowns. Note that all these formulations are "equivalent" , in the sense that one 
can go from one set of unknowns to the other through elementary linear operations (see e.g. ^ ITol/ for 
Maxwell's equations). 

But contrarily to the split PMLs, there is no way of designing unsplit PMLs ( which preserve the original 
operator) for a general first-order hyperbolic system. However, it is straightforward to see that in the 
particular case of a PML corner, with Ux = cTy — cr, the complex coordinate stretching leads to equation 
(O, i.e. to the same equation as the one obtained with the splitting. In that case actually there is neither 
splitting anymore nor additional unknowns. 

3 The model problem: unsplit PMLs for the scalar wave equation 

In the rest of the paper, we focus on a model problem, the two-dimensional acoustic wave equation, written 
as a first-order pressure-velocity system, for which we illustrate the construction of the Cartesian unsplit 
PMLs at the comer and derive some energy estimates. The governing equations of the original first-order 
hyperbolic system in terms of pressure and velocity fields are given by 

-dtP ^ div y, (8) 

p9ty = gradP, (9) 

where P is the pressure, V — (T4, Vy) is the velocity, /i and p are the bulk modulus and the mass density 
respectively, which are assumed to be positive bounded functions. We denote by c = \f[>Jp the acoustic 
sound velocity. Additionally, we should include the initial data for P and V and boundary conditions for 
the pressure field in system (|8]l-(|9), but in order to simplify the presentation, in the rest of the paper they 
will be systematically omitted. 

We introduce the Cartesian unsplit PMLs at the comer following the Zhao-Cangellaris' formulation 
isS^l (recall that the solution of the split Berenger's PML can be deduced by simple linear combinations 
from the solution of the unsplit formulation and conversely, see 1,10,1 ): 
Find (P, P*, y, y*) such that 



-9tP* = divy*, (10) 

p(at + (T,/)y, = a,p, (11) 

p{dt^Oyl)Vy=dyP, (12) 

pdtV:^p{dt^GyI)Vx, (13) 

pdtv; ^ p{dt + oj)Vy, (14) 

-dlP* = -{dt+ a J) [dt + Oyl) P, (15) 
p p 



satisfying the adequate initial conditions, where V* = {V*,Vy) and / is the identity operator In the 
following, we introduce the standard notation, for any arbitrary positive bounded function i/: 

Using the same arguments showed for the Maxwell's equations in ifioll . straightforward computations 
lead to the following result: 

Theorem 2. If <7x ond Uy are constant (but eventually different), the energy 

£2it) = l{\\d?P%^ + \\dfV% + \\axd,V:\\l + \\<yyd,Vy%} , 

satisfies 

j^£2it) = -2 [cTx \\d^V:\\l + ay p^V^X) ^ 
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Finally, the result showed in Ol ill can also be obviously extended here: 



Theorem 3. If cr^ — cry = a, where a is a positive constant, we have the following identity: 



2dt 

where 



F{t) = I P{s)ds. 
Jo 

This implies in particular the decrease of the energy of order 0: 

£o{t) = \{\\P\\l,, + \\vX}- 
4 Semi-discretisation in space using a mixed spectral element method 

The system ([TOll-lfTSll is approximated in space with Qr — Q^^^^ mixed spectral elements based on hexaedral 



meshes, described in 014 



2011 . For the sake of simplicity in its description, we first introduce the spatial 



discretization for the original acoustic model before applying it to the PML problem. 

4.1 Spatial approximation in the fluid domain 

In order to describe briefly the mixed finite elements, we first consider the approximation of the equations 
set in the fluid domain Q, with Dirichlet boundary conditions for the pressure field. The variational formu- 
lation of (O-© is then: 

FindP e H^in), V e such that 

^(y,t/,)p - (gradP,t/,)i2, VtA e {L^mf. 

We introduce T — U^J[Ki a partition of Q with Ne quadrilateral elements, K — [0, 1]^ the unit 
element, and the conform mappings Fi such that Fi{K) — Ki, Vi = 1, . . . , Ne- We set DFi the Jacobian 
matrix of Fi and its Jacobian Ji — det DFi. We finally define the approximation spaces: 

Vla = ['f^C\n), >fi\K,oF,eQr{K), if = on on}, 
v;, - e {L'mf, \J^DF-^xp^^^ o F., e {Qr{K)f] , 

with Qr [K] the set of polynomials whose degree is less or equal to r in each spatial variable. For both 
spaces, the interpolation points coincide with the Gauss-Lobatto quadrature points. The semi-approximate 
problem is then: 

Find Pfi G Pq/j, V G V,^ such that 



dt 
d 



{Vh,'iph)p = {gr ad Ph,xph)L^, yiPh e K- 



Let {(fij : 1 < j < N-p} the finite element basis of Vq/^ and {i/j^^ : 1 < m < N^} the basis func- 
tions of (see ifTil Eofl for a more detailed definition). These basis are associated to the interpolation 
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points that coincide with the quadrature points of the Gauss-Lobatto quadrature formula which is exact for 
polynomials of degree 2r — 1. We introduce the discrete mass and stiffness matrices 

{Rx)jm = (5:r<Pj,'0m)i2 , 1 < j < Nv ■, 1 < m < iVy, 

{Ry)rm = , V„J i2 , 1 < j < A^P, 1 < TO < A^v, 

B,nl = {ip,n,ipi)p, l<m,l<Nv. 

where all the integrals are computed by using the Gauss-Lobatto quadrature formula. The main advantages 
of this method are: (a) it provides diagonal or block diagonal mass matrices (mass lumping) ; (b) the 
stiffness matrices are independent of the mesh and of the physical properties of the fluid medium (gain of 
storage). 

The semi-discretized scheme can then be written in the matrix form: 

j^MPh + R,y^,h + RyVy.h = 0, (16) 
j^BV.^h = RlPh, (17) 
j^BVy,h = R;Ph, (18) 
where we identify the notations for the unknowns and their basis coordinates. 



4.2 Spatial approximation in the PML corner 

Coming back to the semi-discretization of the PML corner problem, we introduce some new matrices, 
defined for any positive function 

The semi-discretized scheme can then be written in the matricial form as: 

^^MP;: + R,V*^ + RyV*^,^0, (19) 
j^B + B''-^V^j, = R*Pi„ (20) 

j^B + B-y^VyM^R;Ph, (21) 



■^MP^ = i —M + -W^+^y + W-^y ) Ph, (24) 



5 Two alternatives for the time discretization 

We introduce At the time step and (P^)", (/'/t)", {VhY+'^, {V*hY^'^) denote the spatial vector fields 
associated to the degrees of freedom of each unknown in the fully discrete problem at time n/S.t and 
(n + 1/2) At, respectively. Since we work with the unsplit PML formulation, one notices that there is 
the second-order equation in time (l24l i to be discretized in the PML comer This leads to two different 
second-order centered finite difference approximations in time. 
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Notation. In the sequel, || • || and (• | •) denote respectively the Euchdian norm and its associated inner 
product. For any positive matrix S, we introduce the notation: 

{U \V)s = {SU \V) ■ \\V\\l = {SV\V) 
5.1 Time Discretization in the fluid domain 

A centered second order finite difference scheme is used for the time discretization of matrix system (fT6l l- 



At 



^ — - ^x^h ) 



B = RyPn , 

completed with the adequate initial conditions. We recall that this scheme is stable under the CFL stability 
condition (e.g. lfll|20|l ): 

{M-^Ku\u) 4 At2 . . , , . 

sup — — < -— r M K IS positive definite, 

where K is the stiffness matrix defined as K = RB^^R*, R is the N-p x 2N\> matrix defined as i? = 
[Rx Ry) and B is the 2A^v x 2A^v block diagonal matrix, each block of size N\> x Ny being equal to B. 
Classically, on a regular grid, this condition is expressed as 

< 1, (25) 

where C is a constant independent of h and At. 

In a homogeneous fluid and using a regular mesh, this condition reduces to 

where d is the dimension of the space and r is the degree of the polynomials. The constant cfld^r can be 
related to the CFL number in dimension one (see Ii20il ) by 

cfld,r = (26) 

In particular for r = 1, we have cfli^i — 1 and in dimension d we recover the classical CFL condition: 

At 1 

which is natural since the Qi — "^'^ discretization on a regular grid is equivalent to the standard second- 
order finite difference scheme. 

5.2 Time Discretization in the PML corner: scheme A 

In this section we consider the scheme which was originally used in [i20i] for the time discretization of the 
governing equations in the PML corners. Let us first introduce some useful notations: 

{DAtU)'' = - u''-^/^)/At ; {lAtU)'' = (f/'^+i/s + t/fc-i/2)/2. 
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The time approximation for (fT9^-(l24l provided by scheme A is written as: 

B{DAtv:^hr = B{DAtv^Mr + s'^^i/a^kjO", 

B{DAtV*hr - BiDAtVy^r + B'^^lAtVy^hT, 

The initial motivation of this study comes from the instabilities observed in the numerical simulations 
for some values of the PML parameters (see Section|6]l. Their origin is located in the PML corner domains 
(both for variable and constant ctq,, a = a;, y), and we also observed that these instabilities were removed 
when decreasing either At or the values of ctq. 

Since the question of stability for variable ctq is an open question even for the continuous PML models, 
we will assume for the analysis of the schemes that a a are positive constants for a = x,y in the corner 
domains. Since ax is constant, we have _B°^^ — axB and we can define the discrete centered operator 
(D'^fU)'' — [DaiU)^ + (Jx{lAtU)'' to approximate dt + (TxI- With these notations and assumptions, the 
scheme can be rewritten as: 



MpAtP,:)"+^ + R.{VlhT+-- + Ry{VluT+-^ - 0, 
B[DZVx,hr = KPh, 

B{D2tVy.hr = RlPJi. 

B{DAtVluT = B[D''^,Vx,Hr, 
B{DAtV*,r = BiDZVy^r, 

MiDl,P;:r = M {{Dl^PhT + [o. + ay){DAtlAtPky 



CTxCTyPh 



(27) 
(28) 
(29) 
(30) 
(31) 
(32) 



We can show the stability result: 



Theorem 4. Assume that ax and ay are constant and equal to a. The discrete scheme A is stable if the 
matrix 



M-—K^ 

is positive definite, where K„ = K + a^M and K = RB^^R*. 

Proof. Since the absorbing functions are equal, it is possible to rewrite the scheme as 

M{DAtPIX+'^ =RAv:jT^'^+Ry{v;jX^K 

B[DAtVi,,r = RiPi:, 

B{DAtVy*jr = R*yPH, 

M{Dl,P*r = M {{Dl.PnT + MDAtlAtPhT 



(33) 



<y^Pl') 



(34) 
(35) 
(36) 
(37) 



Applying D^t to ( l34l l. it is then easy to eliminate V*f-^ by using ( |35] ). 
only in terms of Ph'. 



and dJTb . and rewrite a scheme 



M {{Dl.Phr + MDAtlAtPhT + ^^Ph) + RB-^R*Pi: - 0. 
Introducing the stiffness matrix = RB^^R* + a'^M, we see that the following discrete energy 



£ 



A,h 



pri+1 pn 



At 



{k.Ph\pj: 



M 
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satisfies the identity 




E 



ill 

■A.h 



1<T\\{DAtlAtPhT\\l 



At 



M ' 



and that f ^ ^ > if matrix ( l33T l is positive definite. 



□ 



Let us remark that the sufficient condition stated in the above Theorem for the stability of the numerical 
scheme depends on the value of the absorbing function a, which means that the CFL condition is not the 
same in the PML comer than in the fluid domain. More precisely, if C is the constant appearing in the CFL 
condition in the fluid domain (see (IZSTi), then the stability condition in the PML corner can be expressed as 



which is more restrictive than the one in the fluid domain. 

Application: homogeneous medium and approximation with second order finite differences in time 
and in space. We consider here an homogeneous medium and we use the lowest order mixed spectral 
elements for the approximation, i.e. r = 1, on a regular grid. The scheme can then be written as a finite 
difference scheme in space and time. 



{{DAtfPh)7j = {{{DAtfPhT + K + 'Jv){DAtlAtPhr + ^.^yPh\ 



where the subscripts ij denote the degrees of freedom at nodes of coordinates {iAx, jAy) as usual in the 
finite difference discretizations in the spatial variables. 

The scheme in the fluid corresponds to the classical second order finite difference scheme. In the two- 
dimensional case, its CFL stability condition in the fluid domain corresponds to take C = \/2c in (IZST i and 
so, we obtain 




(38) 




(39) 



In the PML comer, this condition becomes 



cAt < 



h 



1 




(40) 
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5.3 Time discretization in the PML corner: scheme B 

In this section we propose a new scheme which corresponds to another time discretization of the last second 
order differential equation (l24l : 

MpAtP,*)"+^ + + Ry{VluT+'^ = 0, 

B{DAtVy,hr + B"y{iAtVy,hr - R*yPh, 

BiDAtVlhT - BiDAtV^,hr + B^y{lAtV,,hr, 
B{DAtV*hr - B{DAtVy,hr + B'^^ilAtVy.hT, 

For constant values of the damping parameters, this scheme becomes: 

M{DAtP*hT+-- + R.{VlHr+-- + Ry{Vl,y^+-^ = 0, (41) 

S(Z)^==,14,/.)" = RlPJi, (42) 

B^DZVyMT = R;Ph, (43) 

B{DAtV*^^r = B{Dl\V,,^r, (44) 

B{DAtV*^,r = B{DZVy,Hr, (45) 

M{Dl,P*^r = Af (D^==,DX'iP,)", (46) 

In this scheme, the operator {dt + (7xl){dt + i^yP) appearing in (l24l l has been approximated by -Djj^^At' 
which is in some sense more natural than the discretization used in ( |32] |. 
We can show the stability result: 

Theorem 5. We assume that Ux and Uy are constants (not necessarily equal). The discrete scheme B is 
stable if the matrix 

B - —R*M-^R (47) 
4 

is positive definite. 

Proof. We apply D'^Jai to (EB and multiply with D^^P^: 

{DAtlAtMDl^P;, I Dl^P^) + {DljAtVl I R*Dl^Pt) = 0. 
(I) (11) 
We develop each term (I) and (II) and separately: 

(I) = ^ {{DitP^uT^' {Di,p;x-' I iDi,p*,r)^^^, 

and, since <Tx and ay are constants, using ( |46] | and then (|42] |- (|43T |. we obtain 
(II) ^ iDljAtVl\DZDl\R*P,,) 

= (DljAtVl, I D'^^,DZBDZVx,h) + {Dl,lAtVl^ I DZDl^BDZVy,^) 
= iDl,lAtV*n I {DZfDAtV,y)B + [Dl.lAtV*^ I (i^I'-J'i^AtKrjB. 
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The last two terms in (II) can be rewritten as a difference of norms. For instance, the first one leads to 

{Dl,I^tVl^\{DZfD^tV:,,,)B = {Dl,I^tVln\DltV:,,,)B+2a, DIJ^^V*^ [ 

+'jl{Di,i^tVi^\D^al,vi,;)B 



1 

+ 2(7: 

2At 



Hence, if the discrete energy is defined by 



2 












B 





(48) 



then we have the identity 



At 



^2a, ||(i^L/AtV;*,)"||; - 2ay \\{DljAtVy\r\\], < 



In order to point out the stability condition we finally rewrite the discrete energy (l48l l as 



4 1 



2 n*\n+i 



{lAtDAtV*j,r+^^ 

Now, applying DAt to (HTt . it holds 



2 At^ 

M 4 



2 

M 



ilAtDAtV*H) 



* \n+4 



and so 



'-'B,h ~ 



{lAtDi^P*^ 



2 D*\ri+i 



2 At2 , , 

+ {{B - —R*M-'R)iDl,Vir+^ I {Dl,VXr+^) 
M 4 



iiAtDAtv\r+^ 



{lAtDAtv:,) 



Finally, since the matrix (l47l i is assumed positive definite, we obtain £„ , ^ > and conclude the stability 



of the scheme. 



□ 



Corollary 1. // the discrete scheme B is used, then the CFL stability condition Ii47\l holds in the PML 
corner and in the fluid domain. 

Proof. We recall that the CFL condition in the fluid domain is given by requiring M~^RB~^R* is def- 
inite positive whereas the stability condition in the PML comer is obtained by assuming B-^R*M-^R 
is definite positive. 

If we denote by {Xj,Vj), j = 1, . . . , N-p the positive eigenvalues and the eigenvectors of problem 



RB-^R*Vj = XjMv-i, ^ 0, 



(49) 
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the CFL condition in the fluid domain can be expressed as maxj \j < ^p. Analogously, if we denote by 
{lij,Wj), j — 1, . . . , 2N\> the eigenvalues and the eigenvectors of problem 



the CFL condition in the PML corner can be expressed as maxj fij < Then, it is easy to show that 
Xj is a nonzero eigenvalue of i49[ if and only if it is also an eigenvalue of dSOb . In fact, let us assume 
that {Xj, Vj) is an eigenvalue and eigenvector associated to ( |49] |. If we multiply ( l49b by B~^R*M~^ and 
define Wj — B^^R*Vj, we have B^^R* M^^Rwj = XjWj (obviously Wj ^ since Xj and Vj are not 
null). Hence, it is clear that Xj is also an eigenvalue of problem ( l50l l. Reciprocally, we can use analogous 
arguments to show the equivalence between both eigenvalue problems. Hence, maxj fij — maxj Xj and 
consequently the two CFL conditions coincide. □ 



We consider the free propagation of waves generated by a compact supported initial condition (Ricker 
impulse, e.g. ll20ll ') for the pressure field centered at the point source (17, 17), a central frequency equal 
to 1 and a spectral ratio equal to 0.5. The computational domain is [—2,20]^ and the thickness of the 
Cartesian PML is 2 (so that the physical domain is [0, 18]^). The physical parameters are p = 1, ^ = 1. 

The aim of this section is to illustrate numerically the sufficient conditions stated in Theorems |4] 
and |5] with two different spatial discretizations. We consider a spectral finite element method based on 
Qr— Lagrange rectangular piecewise continuous elements for the pressure fields P and P* and piecewise 
discontinuous elements for the velocity fields V and V* on uniform grids. The numerical experiments are 
performed with r — 1 and r = 5. 

Remind that the theoretical CFL condition has been established in the case of a constant absorbing func- 
tion. In the following experiments, we first try to recover these results numerically for constant damping 
functions: 



Then, we consider a quadratic damping function, i.e. we replace the constant value a in ( BTI ) by the 
quadratic profile which is continuous at the inner boundary of the PML and whose upper bound is denoted 
by a*. 

Second-order scheme and constant absorbing function. For r = 1, the positive definite condition ( l47l i 
for scheme B yields to the standard CFL condition ( |39] l, which coincides with the original CFL condition 
for the wave equation and then is independent of the values of and a-y. For scheme A, the CFL condition 
is expected to depend on the value of the absorbing function, as shown in (|40] |. 

In Figures[T]and|2]the continuous lines illustrate the theoretical CFL condition for the discrete schemes 
A and B, for different values of the constant absorbing function (a = 1, 10, 25) and therefore define the 
boundary of the stability region. In both Figures and through the rest of this section, the markers (see the 
circles, triangles and squares in the plots) are located at points {h, At) corresponding to the largest value 
of At for which the schemes have been checked numerically stable in practice. 

In both cases (schemes A and B), the marks lie closed to the boundary of the stability region, so the 
numerical results confirm the predicted CFL condition. For cr = 1, the curve for scheme A is almost the 
same as the one for scheme B, which means that the CFL condition is very closed to the one in the physical 
domain, but this setup corresponds to a very low damping which in practice would require a very large 
PML thickness. As soon as the damping factor is large enough, one can see that the CFL condition is much 
more restrictive than the standard one. 

Higher-order scheme and constant absorbing function. It has been shown in that, for r = 5, 
cfli.5 = 0.1010. According to ( |26] |. we thus have in 2D: C/Z2.5 — 0.1010/a/2. Therefore, since in our 
numerical test c = 1, the CFL condition for the scheme A in the corner PML domain with a constant 



R*M-^Rwj = fijBw. 



w, ^ 0, 



(50) 



6 Numerical illustration 




(51) 
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Figure 1: CFL condition for scheme A with a Figure 2: CFL condition for scheme B with a 
Qi — Q'l^^'^ discretization and a constant absorbing Qi— Qi'^'^ discretization and a constant absorbing 
function. function. 



absorbing function tr, is given by (in with C = ^2/0.1010. In Fi gures |3] and |4] the results are similar 
to the one obtained with r = 1, again confirming the theoretical CFL conditions. In particular the CFL of 
scheme B coincides again with the standard one. However one notices that the CFL condition of scheme 
A is less restrictive than the one obtained with the lower order scheme. 




0.2 0.4 0.6 0.8 1 

h 



Figure 3: CFL condition for scheme A with a 
Qs ^ Q't^'^ discretization and a constant absorbing 
function. 




0.2 0.4 0.6 0.8 1 

h 



Figure 4: CFL condition for scheme B with a 
Q5 ~ Qt""^ discretization and a constant absorbing 
function. 



Quadratic absorbing function. In this paragraph, we consider a continuous quadratic absorbing func- 
tion whose upper bound is given by a* . Notice that the theoretical CFL conditions have been obtained 
for a constant damping, they are therefore not valid anymore in this case. However we compare in this 
paragraph : (i) the stability regions obtained with the theoretical CFL corresponding to a constant damping 
equal to a* and (ii) the numerical stability region for the quadratic profile. 

Figure |6] shows that the CFL condition of scheme B, obtained with Qi — Qf'^*^ finite elements, is still 
independent of a* and thus coincides with the standard CFL (|39] |. These numerical results have been 
checked also for Qs — Q'^^^'^ finite elements, but since they are analogous to those shown in Figure |6] they 
have not been included in the plots. This stability behavior, which was not guaranteed by the theoretical 
results, allows us to conjecture that the CFL of scheme B always remains the standard one whatever the 
damping profile is. 
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Figure |5] (resp. |7} corresponds to numerical experiments performed with scheme A and r = 1 (resp. 
r ~ 5). In both cases the numerical stability regions still depend on the absorbing profile, but this time the 
markers are no longer in the interior of the theoretical stability regions. This could be expected, since the 
theoretical CFL condition has been obtained with the maximum value of the damping profile and therefore 
is more restrictive than the actual CFL condition of the scheme. As in the constant case, the numerical CFL 
condition gets closer to the standard one, when higher-order elements are used. In practice however, since 
we cannot predict the CFL condition for a variable damping, scheme A is not very convenient to use and 
its stability condition is always more restrictive than the one of scheme B. 




Figure 5: CFL condition for scheme A with a Figure 6: CFL condition for scheme B with a 
Qi ~ Qi""^ discretization and a quadratic absorb- Qi — Qf^^'' discretization and a quadratic absorb- 
ing function. ing function. 




A comparison between scheme A and scheme B. Finally, to illustrate the qualitative difference between 
schemes A and B, we use a numerical example where scheme B is stable and scheme A is unstable. To this 
purpose, we present some snapshots of two numerical simulations at different time steps. We have used a 
Qi ^ Qi""^ discretization and a uniform grid with spatial sizes Ax — Ay = 0.5, time step At = 0.2. For 
the construction of the PML, we have used a constant absorbing functions with a — 25. In this case, since 
At < h/ = 0.354, the standard CFL holds in the physical domain and the PML corner domain for the 
scheme B. 
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t=1 t=20 t=30 




Figure 8; Snapshots at time steps n = 1, 20 and 30 for the scheme A. 



As expected, whereas the scheme B remains stable (see Figure |9l), an instability arises by using the 
scheme A (see Figure |8]l. Actually, the more restrictive stability condition for the corner PML, which 
should be used for the scheme A, is here Ai < 0.078. Figure|8]shows that the instability starts at the comer 
PML which is the closest one to the compact support of the initial condition. Moreover this instability 
arises as soon as the wave penetrates the corner PML domain and corresponds to an exponential blow up 
of the numerical solution. 




Figure 9: Snapshots at time steps t — 1, 20 and 100 for the scheme B. 
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Conclusion 

We have emphasized that the PML models in the Cartesian comer domains are always stable and dissipa- 
tive on the continuous level. However, when using the unsplit PMLs, one has to be careful on the time 
discretization of the equation governing the additional unknown stated in the corner domain. In this work 
we have shown that a bad choice for this discretization leads to a scheme which satisfies a restrictive sta- 
bility condition depending on the absorbing function. The numerical solution obtained with this scheme 
blows up exponentially in the comer when one chooses the maximum time step allowed by the stabihty 
condition of the physical domain. This instabiUty coming from the comer is not due to a lack of stability of 
the continuous model and it can be avoided with a right choice of the discretization. This new scheme (la- 
beled as "B") leads to a CFL condition independent of the absorbing function and identical to the stability 
condition of the discrete scheme in the physical domain. 
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